########################################################################
#    Create model script used by BUGS for fitting the 
#    dynamic mixed linear model (dlmm)
# 
#    Author: Wayne Zhang, June 2012 
#            actuary_zhang@hotmail.com
########################################################################

#' Dynamic Linear Mixed Models via Bugs
#'
#' This function fits dynamic linear mixed models in Bugs. 
#' 
#' @param formula  a two-sided linear formula object describing 
#'  the model structure, with the response on the left of a ~ operator 
#'  and the terms, separated by + operators, on the right. The syntax is 
#'  similar to that used by \code{\link[lme4]{lmer}}, where the vertical bar 
#'  character "|" separates an expression for a model matrix and a 
#'  grouping factor for random effects. In addition, the double vertical
#'  bar "||" is used to separate an expression for a model matrix and a 
#'  grouping factor for dynamic effects. Either fixed effects or random 
#'  effects, or both can be missing. However, dynamic effects must be 
#'  present. 
#' @param family  a string in \code{c("gaussian", "bernoulli", "binomial", "poisson")}
#'  that specifies the distribution of the observation.
#' @param data an optional data frame containing the variables named in 
#'  \code{formula}. By default the variables are taken from the environment 
#'  from which the function is called.
#' @param inits  an optional list of initial values to be used for each chain. 
#'  It must be of length \code{n.chains}. Each element is a named list 
#'  with names correponding to the random nodes used in Bugs. If not supplied, 
#'  the function will generate initial values automatically.
#' @param parameters.to.save an optional character vector of the names of the 
#'  parameters to save. 
#' @param weights an optional vector of weights to be used.
#' @param offset this can be used to specify an a priori known component to be 
#'  included in the linear predictor.
#' @param subset an optional vector specifying a subset of observations to be 
#'  used in the creating the model frame.
#' @param na.action a function which indicates what should happen when the data 
#'  contain \code{NA}s. The default is \code{NULL}, where no action is taken. This
#'  is important when the \code{NA}s are used in the response to indicate cells to
#'  be predicted by Bugs, or to flag cells that are truncated or censored. 
#' @param contrasts an optional list. See \code{\link[stats]{lm}}. 
#' @param n.chains number of Markov chains (default: 3). 
#' @param model.file a string that specifies the name of the file storing the 
#'  model script written by the function. This file is further used by 
#'  \code{\link[R2WinBUGS]{bugs}}. The default name is  \code{"model.bug"}. 
#' @param do.fit a logical scalar. When \code{FALSE} the model is not fit but instead 
#'  a structure with the data, initial values and parameter names is returned, which
#'  can be modified for special model forms and passed onto \code{bugs}.
#' @param working.directory sets working directory during execution of this function. 
#'  See \code{\link[R2WinBUGS]{bugs}}. 
#' @param save.ranef a logical scalar indicating whether random effects will be saved.
#'  The default is \code{FALSE}. This is helpful when there is a large number of 
#'  random effects that require considerable amount of memory.
#' @param ... arguments passed to \code{bugs}.
#'
#'
#' @keywords models
#' @export  
#' @include dlmm_bugs_model.R
#' @include dlmm_bugs_data.R
#'   


dlmm_bugs <- function(formula, family = c("gaussian", "bernoulli", "binomial", "poisson"), 
                 data, inits = NULL, parameters.to.save = NULL, 
                 weights, offset, subset, na.action = NULL, contrasts = NULL, 
                 n.chains = 3, model.file = "model.bug", do.fit = TRUE, 
                 working.directory = NULL, save.ranef = FALSE, ...) {
  
  # call
  call <- mc <- match.call()  
  # default na.action in generating model frames
  mc <- as.call(c(as.list(mc), list(na.action = na.action)))
  family <- match.arg(family)
  if (missing(data)) 
    data <- environment(formula)   
  
  # generate (overwrite) data
  data <- bugs_data(mc, formula, contrasts) 
  nc_nlev <- get_ranef_str(mc, formula)
  nc <- nc_nlev$nc   # this could be NULL: no random effects
  nlev <- nc_nlev$nlev
  
  # the dims slot 
  n.obs <- length(data$Y)
  n.beta <- ifelse(is.null(data$X), 0, NCOL(data$X))
  dims <- as.integer(c(n.obs, n.beta,  NCOL(data$F), 
                       length(unique(data$time)), length(nlev),
                       ifelse(is.null(data$wts), 0, 1), 
                       ifelse(is.null(data$off), 0, 1), 
                       family == "gaussian", save.ranef))  
  names(dims) <- c("n.obs", "n.beta", "n.theta", "n.t", "n.term",
                   "n.wts", "n.off", "is.gauss", "save.b")  
  
  # write model script
  if (!is.null(working.directory)) {
    working.directory <- path.expand(working.directory)
  } else {
    working.directory <- tempdir()
  }
  savedWD <- getwd()
  setwd(working.directory)
  on.exit(setwd(savedWD))
  bugs_model(family, dims, nc, nlev, model.file)
    
  # generate initial values if necessary 
  if (is.null(inits))
    inits <- bugs_inits(n.chains, dims, nc, nlev) 
  
  # parameters to be recorded
  if (is.null(parameters.to.save))
    parameters.to.save <- bugs_params(dims)
  
  # output results 
  sims <- list()
  class(sims) <- "bugs"
  out <- new("dlmm", call = call, data = data, inits = inits, 
                parameters.to.save = parameters.to.save, 
                dims = dims, sims = sims)
  if (do.fit)
    out@sims <- bugs(data, inits, parameters.to.save, model.file, 
            n.chains = n.chains, working.directory = working.directory, ...)
  return(out)
}

